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Abstract 

A  two-dimensional  steady-state  model  for  a  phosphoric  acid  fuel  cell  (PAFC)  is  developed.  While  most  of  the  published  literature  deals 
with  one-dimensional  models  for  PAFC,  a  two-dimensional  model  is  necessitated  in  cases  where  the  oxygen  concentration  changes 
substantially  in  the  flow  direction  due  to  depletion  of  oxygen  and  back  diffusion  of  product  water.  The  proposed  modeling  strategy  was 
validated  by:  (i)  testing  in  one-dimensional  mode  for  verification  of  the  basic  parameters  through  a  micro  setup  known  as  “unit  cell”,  and  (ii) 
evaluation  of  the  two-dimensional  model  through  an  experimental  setup  of  a  PAFC  stack  with  four  cells.  Further,  the  utility  of  the  model  in 
design  and  humidity  management,  and  parametric  sensitivity  studies  are  presented. 
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1.  Introduction 

The  need  for  highly  efficient  and  low  emission  energy 
conversion  devices  has  attracted  attention  towards  fuel  cells 
world  over.  A  fuel  cell  is  an  electrochemical  device  which 
converts  chemical  energy  to  electrical  energy  directly.  Out 
of  several  families  of  fuel  cells,  hydrogen-oxygen  fuel  cells 
are  an  important  class  of  fuel  cells  for  which  the  technology 
is  mature  enough  for  practical  usage.  Fuel  cells  have  been 
used  in  the  past  in  space  applications.  Other  possible  appli¬ 
cations  include  small  to  medium  sized  stationary  power 
generation  plants,  vehicle  propulsion,  various  types  of  power 
sources  for  military  use  and  so  on.  Of  the  hydrogen-oxygen 
fuel  cell  systems,  the  most  mature  is  phosphoric  acid 
fuel  cell  (PAFC).  It  operates  at  150-190  °C  and  pressures 
ranging  from  ambient  to  5  atm.  PAFC  systems  primarily 
use  Pt  as  a  catalyst,  both  for  hydrogen  and  oxygen  electro¬ 
des.  Operating  temperature  range  of  PAFC  allows  it  to 
take  up  hydrogen  directly  from  hydrogen  sources  like 
reformer  gases.  CO  (less  than  1%)  present  in  the  reformer 
gases  are  not  adsorbed  on  Pt  sites  owing  to  the  high 
operating  temperature.  The  other  components  used  in  PAFC 
are  mainly  made  of  graphite  and  carbon.  All  these  make 
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PAFC  a  versatile  member  of  the  hydrogen-oxygen  fuel  cell 
family. 

As  a  part  of  the  PAFC  research  program  in  Naval 
Materials  Research  Laboratory  (NMRL),  India,  modeling 
of  PAFC  stacks  of  various  capacities  and  configurations  are 
being  developed.  Numerous  parameters  like  electrode  char¬ 
acteristics,  gas  flow  fields,  acid  concentration,  temperature, 
inlet  partial  pressures  of  the  reactants  and  so  on  have 
significant  effect  on  the  performance  of  a  PAFC  stack.  It 
is,  therefore,  difficult  to  analyze  the  effect  of  each  para¬ 
meter  through  experimentation  in  order  to  design  PAFC 
components  optimally.  A  comprehensive  model  allows  in 
depth  study  of  the  system  at  various  operating  conditions 
with  several  possible  combinations  of  the  component  para¬ 
meters.  Another  advantage  of  a  good  representative  model 
is  easy  development  of  diagnostic  tools  [1]  for  determining 
possible  malfunctions  inside  large  sized,  higher  capacity 
stacks.  Stable  operation  of  medium  to  high  capacity  PAFC 
stacks  require  balancing  of  several  operational  parameters 
like  flow  rates,  humidity  and  temperature.  The  output 
electrical  potential  of  a  PAFC  stack  is  a  function  of  its 
energy  conversion  efficiency.  Thus,  under  variable  load, 
it  is  important  to  maintain  the  potential  at  the  maximum 
possible  value  in  order  to  maximize  energy  conversion 
efficiency.  In  this  paper,  we  develop  a  two-dimensional 
steady-state  model  for  PAFC  to  study  and  understand  the 
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Nomenclature 

aa  effective  area  of  catalyst  site  per  unit  volume 

of  agglomerate  (cm2  cm-3) 
ac  groove  cross-section  (cm2) 

A*  cross-section  area  of  an  agglomerate  in 

FZ-plane  (cm2) 

b  normal  Tafel  slope  (V  per  decade) 

cos  concentration  of  oxygen  on  surface  of 

agglomerate  (gmol  cm-3) 

c  ia  concentration  of  oxygen  inside  agglomerate 

(gmol  cm-3) 

ch3po4  phosphoric  acid  concentration  in  percent  ortho 
phosphoric  acid  equivalent 
Dagg  effective  diffusivity  inside  agglomerate  for 
oxygen  (cm2  s-1) 

Dd  effective  diffusivity  in  diffusion  layer  for 

oxygen-water  (cm2  s-1) 

Dr  effective  diffusivity  in  reaction  layer  for 

oxygen-water  (cm2  s-1) 

Dhc  hydraulic  diameter  of  each  groove  (cm) 

E  open  circuit  potential  (V) 

E°  standard  open  circuit  potential  (V) 

El  oxygen  electrode  (solid  carbon  phase) 

potential  (V) 

E2  local  electrolyte  (liquid)  potential  (V) 

/  friction  factor 

F  Faraday  constant  (96,500  C  g-1  equivalent 

weight) 

H  Henry’s  constant  for  oxygen  solubility  in 

phosphoric  acid  (gmol  cm-3  atm-1) 

4  local  current  density  (A  cm-2) 

i0  exchange  current  density  for  oxygen  reduction 

on  Pt  (A  cm-2) 

ix  current  in  X- direction  (A) 

j(x)  flow  of  cation  in  X-direction  (A  cm-2) 

j(y)  flow  of  cation  in  7-direction  (A  cm-2) 

k  reaction  constant  of  oxygen  reduction  reaction 

lc  groove  length  (cm) 

mi  product  of  molar  flux  of  ith  species  at  gas 

electrode  interphase  in  X-direction  and  width 
of  the  cell  (gmol  cm-1  s-1) 

niio  inlet  molar  flow  rate  of  ith  species  (gmol  s-1) 

n  no.  of  electrons  taking  part  in  the  reaction 

ns  grooves  per  manifold 

Al  no.  of  agglomerate  exposed  per  unit  sqcm  area 
in  FZ-plane 

Pi  partial  pressure  of  the  ith  species  (atm) 

Pto  partial  pressure  of  the  ith  species  at  inlet  (atm) 
Pi  total  pressure  at  groove  inlet 

R  universal  gas  constant  (atm  cm3  gmol-1  k-1) 

g  average  radius  of  an  agglomerate  (cm) 

Ron  oxygen  reduction  rate  (gmol  s-1  cm-3) 
r  radial  dimension  inside  an  agglomerate 

originating  at  the  center  (cm) 


T  system  temperature  in  kelvin  scale 

U  oxygen  reaction  potential  wrt  the  reference 

(V) 

W  width  of  the  cell  in  Z-direction  (cm) 

WAax  maximum  electrical  work  available  (W) 

A ,£*?ey  standard  half  cell  potential  wrt  hydrogen 

electrode  (V) 

Greek  letters 

f]  local  overpotential  (V) 

rj0  overpotential  at  the  bulk-electrolyte  and  elec¬ 

trode  interphase  (V) 

Ke  electrolyte  conductivity  (mho  cm-1) 

pa  agglomerate  density  inside  reaction  layer 

(no.  cm-3) 

W  flow  per  manifold  at  system  temp  and  pressure 

(cm3  s-1) 

p  viscosity  (g  cm  1  s_1) 

Subscripts 

1  oxygen  (ith  species) 

2  water  (ith  species) 

relationship  between  the  performance  of  the  fuel  cell  and 
various  design  options. 

2.  Phosphoric  acid  fuel  cell  (PAFC) — principle 
of  operation 

A  phosphoric  acid  fuel  cell  is  composed  of  two  porous  gas 
diffusion  electrodes  namely  anode  and  cathode  (see  Fig.  1) 
[2]  juxtaposed  against  a  porous  electrolyte  matrix.  The  gas 
diffusion  electrodes  comprise  a  porous  substrate  that  faces 
the  gaseous  feed.  This  substrate  is  a  porous  carbon  paper  or 
cloth.  On  the  other  side  of  this  substrate,  which  faces  the 
electrolyte  (phosphoric  acid),  platinised  fine  carbon  powder 
electrocatalyst  is  roll  coated  with  poly-tetra-fluro-ethylene 
(PTFE)  as  a  binder.  PTFE  also  acts  as  a  hydrophobic  agent  to 
prevent  flooding  of  pores  so  that  reactant  gas  can  diffuse  to  the 
site  easily.  The  two  electrodes  are  separated  by  phosphoric 
acid  (electrolyte;  Fig.  1).  At  the  anode,  hydrogen  ionizes  to 
H+  and  migrates  towards  the  cathode  to  combine  with  oxy¬ 
gen,  forming  water.  The  product  water  then  diffuses  out  to  the 
oxygen  stream  and  comes  out  of  the  system  as  steam.  An  emf 
is  generated  between  the  two  electrodes  through  conversion 
of  reaction  free  energy  to  electricity  and  on  connecting  an 
external  load  electrical  power  can  be  extracted. 

The  reactions  at  anode  and  cathode  are  as  follows  [2]. 
At  anode: 

H2^2H++2e-,  A£?ev=0  (1) 

At  cathode: 

\  02  +  2H+  +  2e-  =  H20,  A£r°ev  =  1 .229  (2) 

02  +  2H+  +  2e“  =  H202,  A£r°ev  =  0.67  (3) 
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Fig.  1.  Phosphoric  acid  fuel  cell — principle  of  operation. 


H202  +  2H+  +  2e"  =  2H20,  A£?ev  =  1 .7  (4) 

H202  ^H20  +  i02  (5) 

Under  reversible  condition  maximum  electrical  work 
available  is  [2]: 

lUmax  —  A G water  formation  (6) 


and  the  ideal  open  circuit  potential  is  given  by  Nernst 
equation  [2]: 


E 


AG 

nF 


=  £°  + 


^hydrogen  Oxygen 

^water 


In  Fig.  2,  an  exploded  view  of  a  practical  single  cell  is 
depicted  [2]  where  syrupy  phosphoric  acid  is  immobilized 


exhaust 


Fig.  2.  Exploded  view  of  PAFC  cell. 
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inside  a  SiC  matrix.  The  high  density  graphite  plates  as 
shown  in  Fig.  2  have  grooves  on  one  side  through  which 
reactant  gases  are  provided  to  the  electrodes  (anode  and 
cathode).  The  porous  substrate  of  the  electrodes  face  the 
grooved  side  of  the  graphite  plates,  whereas  the  catalyst  is 
coated  on  the  other  side.  In  between  the  two  electrodes, 
phosphoric  acid  is  provided  in  silicon  carbide  matrix. 

The  whole  assembly  is  then  sandwiched  between  two 
insulator  plates  (not  shown  in  the  figure)  and  bolted  tightly. 
The  two  graphite  plates  act  as  the  terminals  for  connecting 
the  electrical  load.  Such  cells  are  stacked  in  series/parallel 
combinations  for  useful  power.  In  large  stacks,  the  graphite 
plates  are  bipolar  in  nature  where  one  side  of  a  plate  provides 
hydrogen  to  one  cell  and  the  other  side  oxygen  to  the 
adjacent  cell  in  series  [3]. 

3.  Background  of  model  development 

Modeling  of  a  phosphoric  acid  fuel  cell  requires  the 
representation  of  the  following  factors  in  the  model:  (i) 
temperature,  (ii)  operating  pressure,  (iii)  humidity,  (iv) 
electrical  load,  and  (v)  phosphoric  acid  management  inside 
PAFC. 

PAFC  electrodes  are  basically  porous  gas  diffusion  elec¬ 
trodes  (GDE),  which  are  difficult  to  characterize  on  a 
microscopic  level.  One  approach  to  this  problem  is  to 
employ  a  macroscopic  model  like  “porous  electrode  the¬ 
ory”  [1,4]  that  accounts  for  the  essential  features  of  a  porous 
electrode  without  going  into  the  exact  geometric  details. 
GDEs  like  that  of  a  PAFC  adds  an  additional  degree  of 
complexity  due  to  presence  of  three  phases,  viz.  gas,  elec¬ 
trolyte  and  catalyst  surface.  This  requires  a  model  with  dual 
porosity,  one  each  for  the  gas  phase  and  the  electrolyte 
phase.  One  approach  is  to  employ  “flooded  agglomerate” 
concept,  first  introduced  by  Giner  and  Hunter  [5].  A  mod¬ 
ified  version  of  this  concept  was  exploited  for  the  develop¬ 
ment  of  the  model  in  this  work  as  described  later.  The  main 
advantage  of  considering  this  concept  is  its  easy  adaptability 
to  electrode  made  by  various  methods.  It  should  be  recog¬ 
nized  that  the  flooded  agglomerate  model,  when  executed 
correctly,  is  essentially  a  discretized  version  of  “porous 
electrode  theory”  [1]. 

Fig.  2  shows  the  cross-section  of  a  composite  cathode. 
Bulk  gas,  i.e.  oxygen,  flows  in  the  T-direction  through  the 
ribbed  support  over  the  left  side  of  the  electrode.  As  the 
gas  flows,  part  of  the  oxygen  diffuses  through  the  porous 
substrate  of  the  electrode.  The  diffusion  of  oxygen  is  in 
X-direction.  The  gas  after  passing  through  porous  substrate 
(diffusion  layer)  reaches  the  catalyst  layer  also  termed  as 
reaction  layer.  Catalyst  layer  is  a  porous  layer  partly  filled 
with  electrolyte.  The  gas  diffuses  through  the  dry  pores, 
dissolves  in  the  electrolyte  and  then  diffuses  to  the  catalytic 
site.  Oxygen  reduction  takes  place  on  the  site  and  the 
product  water  back  diffuses  through  the  diffusion  layer  to 
escape  in  the  oxygen  stream. 


Most  of  the  published  literature  deals  with  one-dimen¬ 
sional  models  for  PAFC  cathode.  This  is  valid,  provided  the 
throughput  of  the  gas  is  quite  high  inside  the  cell,  so  that 
the  change  in  concentration  of  the  gas  in  T-direction  is 
negligible.  However  for  a  big  cell  with  limited  throughput 
of  gas,  concentration  of  oxygen  changes  substantially  in 
y-direction  owing  to  depletion  of  oxygen  and  back  diffusion 
of  product  water.  This  leads  to  the  necessity  for  developing  a 
two-dimensional  model  for  the  oxygen  electrode,  consider¬ 
ing  no  changes  in  Z-direction.  This  condition  is  valid  for 
most  of  the  cell  as  the  gas  feeding  is  done  by  parallel  grooves 
along  Z-direction. 

It  is  generally  accepted  that  the  major  difficulty  in  fuel  cell 
development  lies  with  the  kinetics  of  the  oxygen  electrode 
[2,6].  The  hydrogen  electrode  kinetics  is  relatively  much 
faster  than  that  of  the  oxygen  electrode.  Most  of  the 
researchers  modeling  fuel  cells  have  considered  primarily 
two  aspects,  viz.  diffusion  of  gases  to  electrode  surface  and 
kinetics  of  the  oxygen  electrode.  To  the  best  knowledge  of 
the  authors,  two-dimensional  models  for  PAFC  cathode  in 
the  literature  have  been  minimal.  But  many  two-dimensional 
models  have  been  proposed  for  other  types  of  fuel  cells. 
Typically,  two-dimensional  models  are  developed  to  study 
and  analyse  various  issues  concerning  the  fuel  cell,  e.g. 
electrocatalyst  utilization,  current  density  distribution, 
hydrodynamics,  etc.  Boersma  and  Sammes  [7]  have  ana¬ 
lysed  the  distribution  of  gas  flow  in  internally  manifolded 
solid  oxide  fuel  cells  (SOFC).  Kulikovsky  et  al.  [8]  have 
considered  a  two-dimensional  model  to  study  the  dead  and 
active  reaction  zones  in  the  cathode  compartment  of  a 
polymer  electrolyte  membrane  (PEM)  fuel  cell.  Dutta 
et  al.  [9]  have  used  a  three-dimensional  flow  model  to 
analyse  the  effects  of  membrane  thickness  and  cell  voltage 
on  the  axial  distribution  of  current  density  and  rate  of  water 
transport  in  PEM  fuel  cells. 

A  number  of  two-dimensional  models  for  proton 
exchange  membrane  fuel  cell  (PEMFC)  have  also  been 
published,  the  applications  of  which  have  been  varied — to 
study  water  transport  [10];  influence  of  ohmic  drop  within 
the  active  layer  on  electrode  performance  [11];  calculation 
of  cell  performance,  ohmic  resistance  and  water  profile  in 
the  membrane,  current  distribution  and  variation  of  tem¬ 
perature  along  the  gas  channels  [12];  hydrodynamics  of  the 
gases  and  multicomponent  transport  in  the  cathode  [13];  a 
computational  fluid  dynamics  (CFD)  based  model  to 
account  simultaneously  for  electrochemical  kinetics,  current 
distribution,  hydrodynamics,  multicomponent  transport  and 
hydrogen  dilution  effects  [14];  two-dimensional,  two-phase, 
multicomponent  transient  model  [15]. 

Some  of  the  above-mentioned  issues  are  also  studied  for 
direct  methanol  fuel  cells  (DMFC)  [16-18]  and  molten 
carbonate  fuel  cells  (MCFC)  [19-21].  To  undertake  a  com¬ 
prehensive  study  and  a  detailed  analysis  of  PAFCs,  as  a  first 
step  in  the  process  engineering,  we  present  a  two-dimen¬ 
sional  steady- state  model  in  this  work.  The  model  has  been 
used  to  study  the  effects  of  various  process  parameters  on  the 
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performance  of  PAFC  and  gain  insight  into  the  different 
phenomena  occuring  in  PAFCs  to  help  in  carrying  out 
further  advanced  studies  like  control  and  diagnostics. 


4.  Development  of  steady  state  model  for  PAFC 
cathode  (oxygen  electrode) 


Diffusion  of  Oxygen  and  back  diffusion  of  water  vapor  in 
the  porous  substrate,  where  no  reaction  takes  place,  can  be 
defined  by  the  following  governing  equation: 


■DA  (fpA  (da  (&Pi\ 

RTJ  \dx 2  )  \RTj  \dy2  ) 


f  =  1,2 


Modeling  of  the  reaction  layer  can  be  done  by  using  the 
flooded  agglomerate  model  [1,4]  coupled  with  gas  diffusion 
equations. 

Fig.  3  depicts  the  concept  of  flooded  electrolyte  model  [1]. 
It  is  considered  that  the  carbon  particles  containing  catalyst 
(Pt)  formed  small  agglomerates.  These  agglomerate  glo¬ 
bules  are  completely  flooded  with  electrolyte  (phosphoric 
acid).  The  inter  agglomerate  space  is  partially  filled  with 
electrolyte  due  to  the  presence  of  PTFE.  Reactants  (oxygen) 
can  diffuse  through  this  inter  agglomerate  space  and  can 
reach  the  agglomerate  surface.  The  reactant  is  consumed  in 
the  agglomerates  and  product  water  diffuses  out  from  the 


agglomerate  surface.  Diffusion  of  gas  in  electrolyte  medium 
is  inside  the  agglomerates  only.  Although  the  size  of  the 
agglomerate  particle  may  be  somewhat  uncertain  and  may 
vary  with  position,  still  this  model  gives  a  good  account  of 
the  processes  involved  and  can  be  adopted  easily  to  different 
types  of  electrodes. 

The  following  assumptions  are  made  in  the  model  devel¬ 
opment  [1,4]: 


•  Spherical  agglomerates. 

•  Agglomerate  size  <C  thickness  of  the  active  reaction  layer. 

•  No  ohmic  loss  inside  carbon  particles. 

Hence,  oxygen  transport  inside  the  spherical  agglomerate 
undergoing  reaction  can  be  written  as  [1]: 


d_ 

dr 


with  the  following  boundary  conditions  [1]: 


^la  —  ^os 


(10) 


(ii) 


Since  all  fuel  cell  operations  are  carried  out  much  away 
from  equilibrium,  that  is  the  overpotential  is  always 


support 


Fig.  3.  Flooded  electrolyte  model. 
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>75  mV,  Ron  in  equation  (9)  can  be  assumed  to  exhibit  Tafel 
kinetics: 


(12) 


where  k  =  i0/  ( nFcos ). 

Eq.  (9)  and  (12)  along  with  the  boundary  conditions, 
Eqs.  (10)  and  (11),  can  be  analytically  solved  to  give  [1,22]: 


c  h 


c 


os 


C^agg/r) 


sinh  (4>(r)) 
sinh  (0(flagg)) 


where 


(13) 


rexp 


(14) 


The  molar  flux  of  oxygen  (Noxy)  at  the  surface  of  an 
agglomerate  can  be  estimated  from  Eq.  (13)  as 


^agg Cos  [</>  (^agg  )  COth ( (j)  (R^gg  )  )  1  ] 


(15) 


Henceforth,  consumption  of  oxygen  per  unit  volume  at 
any  point  inside  the  reaction  layer  can  be  determined  by  the 
product  of  agglomerate  density,  single  agglomerate  surface 
area  and  Noxy.  Therefore  the  governing  equation  for  oxygen 
balance  in  the  reaction  layer  can  be  written  as 


'PA  (&pA  (pA  (fpf\ 

RTJ  \  dx2  J  \RTJ  \dy2  J 


(16) 


where 


r  47rpaZ)agg^aggC0S  [(j)  (flagg)  COth(0  (/?agg ) )  1]  (12) 


0(^agg)  — 


and 


(18) 


Cos 


(19) 


Distribution  of  overpotential  inside  the  catalyst  layer  can 
be  estimated  by  considering  ohmic  law  for  the  cation  flow,  in 
this  case  hydrogen  ions,  in  the  electrolyte  filled  zone  [4]. 
Thus,  change  of  overpotential  in  X-and  y-direction  can  be 
related  to  the  current  flux  as  follows: 


2i=  iff 

dx  kqA*N  1 

dy=  j(y) 

dy  KqA*N1 


(20) 

(21) 


Considering  a  small  control  volume  of  Ax  and 
Ay  (lim  Ax,  Ay  — >  0)  inside  the  catalyst  layer  and  balancing 


the  charge  flow  across  it  yields 


dj(x)  djjy) 
dx  dy 


=  —nFF 


(22) 


Differentiating  Eqs.  (20)  and  (21)  and  substituting  them  in 
Eq.  (22)  gives 

d2f]  d2f]  nFF 

— -  4 - -  = -  (23  ) 

dx2  dy2  KeA*Nl 


Cation  flow,  i.e.  hydrogen  ion  flow  in  the  catalyst  layer 
depends  upon  the  local  electrolyte  conductivity  which  in 
turn  depends  upon  the  local  phosphoric  acid  concentration. 
A  non-linear  curve  was  fitted  to  co-relate  electrical  con¬ 
ductivity  and  phosphoric  acid  concentration  from  the  data 
available  in  the  literature  [23]  as  shown  in  Eq.  (24). 

Ke  =  0.72796  —  6.5068  exp(—07)cH3po4  (24) 


It  is  assumed  that  the  concentration  of  acid  inside  the 
agglomerates  remain  constant  and  the  local  partial  pressure 
of  the  acid  immediately  outside  the  agglomerate  is  in 
equilibrium  with  the  acid  inside  the  agglomerate.  Consider¬ 
ing  this,  the  phosphoric  acid  concentration  and  equilibrium 
vapor  pressure  at  system  temperature  was  co-related  from 
the  data  available  in  the  literature  [24] .  The  best  fitted  curve 
is  shown  in  Eq.  (25). 


Ch3po4 


1 

0.008623  +  0.004187^2 5 


(25) 


Water  transfer  inside  the  catalyst  layer  can  be  depicted  by 
considering  the  fact  that  local  water  generation  due  to 
reaction  diffuses  out  towards  the  bulk  gas  alongwith  eva¬ 
porated  moisture  from  the  bulk  acid  catalyst  layer  inter¬ 
phase.  Thus,  local  acid  concentration  changes  depending 
upon  the  position  from  the  bulk  gas  electrode  interphase  and 
the  local  current  density.  Higher  current  density  means  more 
dilution  of  acid,  however,  proximity  to  bulk  gas  zone  or 
diffusion  layer  implies  drier  environment.  Considering  neg¬ 
ligible  mass  transfer  resistance  from  bulk  to  agglomerate 
surface  due  to  intimate  positioning,  partial  pressure  of  water 
vapor  at  agglomerate  surface  remains  same  as  that  of  the 
bulk  phase  in  the  catalyst  layer.  Therefore,  the  water  trans¬ 
port  phenomena  in  the  catalyst  layer  can  be  modeled  by  the 
following  equation: 


PA  (fpf  |  (da  (fiA 

RTJ  \  dx2  J  \RTj\dy 2  J 


=  2  r 


(26) 


4.1.  Boundary  conditions 

Oxygen  diffusion  rate  through  the  porous  layer  is  deter¬ 
mined  by  the  rate  of  consumption  of  the  gas  in  the  catalyst 
layer  as  well  as  the  concentration  of  the  gas  at  bulk. 
Depending  upon  the  flow  rate,  bulk  oxygen  concentration 
changes  along  the  T-direction  (Fig.  2).  Also  back  diffusion 
of  water  vapor  changes  the  concentration  profile  inside  the 
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porous  substrate.  Taking  into  consideration  all  these  effects, 
the  system  can  be  solved  as  a  boundary  value  problem. 
However,  concentration  profile  along  the  bulk  flow  direction 
is  not  known,  as  it  is  determined  by  the  inlet  flow  rate, 
composition,  pressure  drop  along  the  bulk  flow  channel  and 
the  diffusion  rates  of  the  reactant  and  product.  To  overcome 
this  problem  an  initial  guess  value  of  the  concentration 
profile  in  the  bulk  could  be  given.  After  one  set  of  computa¬ 
tion  is  over,  this  profile  may  be  corrected  from  the  diffusive 
fluxes  at  different  positions  in  the  bulk  flow  direction.  This 
process  may  be  continued  till  the  profile  does  not  change 
further. 

Considering  a  uniform  pressure  profile  as  initial  guess  for 
oxygen  and  water  vapor: 

Pi=Pio,  atx  =  0,  Vy 


P2=P20,  atv  =  0  ,Vy 


After  one  set  of  computation,  dpt/dx  at  v  =  0,Vy  is 
determined.  Considering  a  control  length  of  Ay  on  the 
bulk-porous  layer  interface  and  taking  mass  balance  over 
it  and  defining  mt  as  the  product  of  molar  flux  of  ith 
component  at  the  bulk  gas-electrode  interface  and  width 
of  the  electrode,  W  (refer  Fig.  2): 


nti  =  W 


i  =  1,2 


(27) 


Considering  pressure  drop  along  the  flow  channel  due 
to  hydraulic  resistance  (Fig.  2),  partial  pressure  along  T-at 
X  =  0  can  be  rewritten  as 


Pi\x=0,y  Pt 


fo  rn  dy  +  mi0 

E;=  i  fo  mi dv  +  mj0 


i  =  1,2 


(28) 


In  order  to  find  out  the  local  total  pressure  along  the  bulk 
gas  electrode  interphase  it  may  be  considered  that  the  flow 
inside  the  channel  is  in  laminar  zone.  This  is  generally 
valid  for  the  planar  structure  fuel  cells  upto  a  flow  length 
of  1  m  [25].  It  can  be  also  checked  by  determining  the  actual 
Reynolds  number  for  specific  cases.  Thus,  considering 
Poiseuille  flow: 


A p  =  KcW  (29) 

where  Kc  =  fplc/nsacD^c. 

W  being  local  flowrate  per  cell  at  system  temperature 
and  pressure;  /  is  the  friction  factor  which  depends  on  the 
channel  geometry,  i.e.  /  =  28.4  for  square  cross-section 
[25].  Thus,  after  every  set  of  computation,  the  partial 
pressure  of  water  vapor  and  oxygen  can  be  updated  along 
the  bulk  flow  direction  till  the  profile  remains  unchanged.  To 
provide  the  boundary  conditions  for  water  transport  equation 
inside  the  electrode,  it  can  be  considered  that  at  electrode 
and  bulk  acid  interphase  the  gas  remains  saturated  with  bulk 
acid  concentration. 

At  v  =  B,  Vy: 

P2  =  Plsat  (BC  1 ) 


As  shown  in  Fig.  2,  at  inlet  entire  electrode  cross-section 
is  exposed  to  the  incoming  stream. 

At  y  =  0,  Vv: 

Pi  =  Pio  (BC  2) 

P2  =  P20  (BC  3) 


At  y  =  L,  Vx 

dpi 


dy 

dpi 

dy 


0  (BC4) 


=  0  (BC5) 


Boundary  conditions  4  and  5  (BCs  4  and  5)  are  justified  by 
the  fact  that  at  the  end  of  the  electrode  entire  cross-section  is 
blocked  with  the  sealant  and  no  material  transfer  takes  place 
across  the  cross-section. 

At  x  =  B,  Vy: 

tr=°  <BC6> 


BC  6  is  valid  considering  no  oxygen  transport  beyond  the 
electrode’s  catalyst  layer.  The  reference  potential  is  the 
hydrogen  electrode  which  is  placed  after  the  catalyst  layer 
of  the  oxygen  electrode  separated  by  the  acid  electrolyte. 
Hence,  the  overpotential  measured  in  reference  to  it  can  be 
considered  as  the  overpotential  at  the  end  of  the  oxygen 
electrode  catalyst  layer  (BC  7). 

At  x  =  B,  Vy: 

V  =  Vo  (BC  7) 

At  v  =  A,  Vy: 


drj 

dx 


0  (BC8) 


At  y  =  0,  A  <  v  <  B: 
n  =  '/o  (BC9) 

BC  8  is  considered  since  ion  movement  takes  place  only 
inside  the  electrolyte  filled  catalyst  layer  and  not  in  the  dry 
diffusion  layer. 

At  y  =  L,A  <  x  <  B: 

1^  =  0  (BC  10) 
dy 


4.2.  Model  solving  and  validation 

The  model  developed  was  numerically  solved  using 
Pdease2D.  Pdease2D  is  a  finite  element  software  for  solving 
two-dimensional  linear/non-linear  partial  differential  equa¬ 
tions.  Pdease2D  solves  simultaneous  PDE  using  Galarkin 
finite  element  method.  It  has  adaptive  grid  control  mechan¬ 
ism  built  in.  Since  the  electrode  structure  is  composed  of 
two  layers  it  was  so  programmed  that  the  variables  are 
continuous  at  the  layer  interfaces.  Macsyma  programming 
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language  platform  was  used  for  customizing  and  interfacing 
with  peripheral  software. 

The  set  of  equations  are  solved  with  a  given  overpotential 
at  the  electrolyte  side.  Further,  integrating  the  molar  flux  of 
oxygen  in  X-direction  along  the  bulk  gas-electrode  interface 
and  multiplying  by  width  of  the  cell,  W,  total  oxygen  inflow 
can  be  found.  Once  the  molar  flux  of  oxygen  has  been 
determined,  total  current  generated  can  be  calculated  as 

pi  c 

gross  current  =  nFW  /  m\  dy  (30) 

Jo 

Thus,  the  steady  state  overall  current  generated  against  a 
given  overpotential  and  oxygen  flow  rate  with  inlet  humidity 
can  be  determined. 

The  model  developed  was  first  tried  in  one-dimensional 
mode  to  see  the  validity  of  the  basic  parameters.  In  order  to 
test  this  a  micro  setup  known  as  “unit  cell”  was  used  where 
the  flowrate  kept  is  such  that  there  is  no  concentration 
gradient  along  the  bulk  flow  interphase.  After  this  further 
tests  were  carried  out  in  experimental  stacks,  details  of 
which  are  provided  in  the  following  section  and  difference 
between  one-dimensional  and  two-dimensional  model  was 
analysed. 

4.2.1.  Experimental  setup  for  model  validation 

To  validate  the  model  in  one-dimension,  i.e.  “X-direction” , 
the  “unit  cell”  setup  is  used.  It  consists  of  a  single  cell  of  one 


cathode  and  one  anode  separated  by  acid  layer.  The  area  of 
the  exposed  electrode  is  very  small  so  that  the  concentration 
of  the  gas  at  any  point  over  the  electrode  surface  remains 
same.  The  electrodes  were  connected  through  a  “PARC” 
potentiostat  galvanostat  with  EG  &  G  developed  electro¬ 
chemical  experimentation  software.  The  sketch  of  the 
experimental  setup  is  shown  in  Fig.  4.  A  hydrogen  electrode 
is  used  as  the  reference  point  for  the  measurement.  This 
arrangement  allows  independent  readings  of  each  electro¬ 
des.  The  detail  design,  being  proprietary  of  NMRL,  is  not 
discussed  here. 

To  validate  the  model  in  two-dimensional  form,  an  aver¬ 
age  of  four  identical  experimental  phosphoric  acid  fuel  cell 
stack  each  consisting  of  four  cells  were  used.  In  effect  an 
overall  average  of  16  cells  separated  in  four  groups  were 
used  for  model  validation.  All  the  components  of  the  stack 
were  developed  in  house  at  NMRL. 

A  conceptual  sketch  of  the  stack  and  the  components  are 
shown  in  Fig.  5.  As  shown  in  the  figure,  two  types  of  groove 
design  were  used,  types  A  and  B.  The  electrodes  were  made 
by  roll  coating  catalyst  powder  on  hydrophobised  porous 
TGP-120  graphite  paper.  The  catalyst  used  was  Pt  loaded  on 
Vulcan-XC72R  carbon.  Teflon  suspension  of  particle  size 
range  *-*  was  used  as  a  binder  and  hydrophobic  agent  to 
prevent  flooding  of  pores  in  the  catalyst  layer.  In  between  the 
two  electrodes,  a  silicon  carbide  matrix  was  incorporated  to 
hold  syrupy  concentrated  phosphoric  acid. 


Fig.  4.  Unit  cell  experimentation  setup. 
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TYPE  —A  grooves 


The  manifold  plates  on  the  back  side  is  similarly  connected  to  the  hydrogen  headers 

Bipolar  graphite  gas  manifold  plate 


Scehematic  of  stack  assembly 

Fig.  5.  Experimental  stack  components  and  assembly  scheme. 


In  a  stack,  all  the  cells  are  electrically  connected  in  series 
and  the  reactant  gases,  viz.  humidified  hydrogen  and  oxygen 
are  provided  through  common  headers  to  each  cell.  Simi¬ 
larly,  all  the  outlets  are  connected  to  another  common  header 
through  which  the  lean  gas  flows  out.  In  order  to  eliminate 
maldistribution  of  gases  in  the  cells  in  a  stack,  fewer  number 
of  cells  (four  cells)  were  assembled  with  extra  large  feeder 
holes.  The  assemblage  was  done  like  a  plate  and  frame  filter 
press  assembly.  Conceptual  flow  pattern  of  gases  inside  the 
stacks  is  shown  in  Fig.  5.  Each  cell  is  characterised  as 

•  effective  electrode  area  =  100  cm2; 

•  groove  type  A  =  straight  grooves,  all  parallel,  gas  travel 
length  =10  cm,  width  =10  cm; 

•  groove  type  B  =  groove  tortuous  and  are  carried  over 
from  inlet  to  outlet,  gas  travel  =  20  cm,  width  =  5  cm. 


Schematic  of  the  complete  experimental  setup  is  shown  in 
Fig.  6.  Pure  hydrogen  and  oxygen  gases  were  used  from  the 
cylinders.  The  hydrogen  was  humidified  with  predetermined 
humidity  level,  in  order  to  prevent  drying  up  of  acid,  while 
the  oxygen  was  used  directly.  In  the  outlet  streams,  humidity 
levels  of  hydrogen  and  oxygen  was  checked  using  dry  and 
wet  bulb  methods.  Emf  of  all  the  cells  were  measured  using 
precision  voltmeter  and  the  circuit  current  was  measured 
using  shunt  based  ammeter.  Experimentation  was  done 
by  keeping  the  overall  stack  voltage  at  a  predetermined 
value.  In  order  to  operate  the  system  at  such  potentiostatic 
mode  a  dc  compensator  was  used  in  series  with  the  stack. 
An  average  of  all  the  four  cells  in  series  is  taken  as  the  cell 
potential. 

Parameters  like  agglomerate  radius,  agglomerate  density 
are  determined  through  scanning  electron  micrographs  of 
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the  electrodes.  Other  parameters  like  “Tafel  slope”,  and 
active  area  per  unit  volume  of  agglomerate  etc,  are  tuned  to 
match  the  model  output  to  some  standard  experimental 
results  from  small  test  PAFC  stacks.  Table  1  lists  the  final 
values  of  the  parameters  used  for  simulating  the  actual 
stacks. 

4.2.2.  Analysis  of  data  for  model  validation 

Experimental  runs  were  taken  on  the  unit  cell  to  check 
the  model’s  one-dimensional  prediction.  To  simulate  one¬ 
dimensional  effect  the  concentration  along  the  bulk  gas 
flow  is  kept  unchanged.  Fig.  7  depicts  experimental  and 
simulated  values  at  various  overpotentials.  The  values  match 
quite  well  and  as  expected  the  current  increases  rapidly  with 
overpotential  as  the  reaction  rate  increases  exponentially 
with  respect  to  the  overpotential,  refer  Eq.  (12). 

Several  experimental  runs  were  taken  under  various  flow 
rates  and  output  potential  settings.  On  the  hydrogen  side, 
flow  rate  was  always  kept  twice  the  oxygen  flow  rate  to 
maintain  adequate  pressure  differential.  Temperature  of  the 
fuel  cell  was  maintained  at  150  °C.  The  oxygen  inlet  was 
maintained  at  0%  humidity  and  the  hydrogen  inlet  was 
saturated  with  moisture  at  60  °C.  The  reason  for  maintaining 
this  particular  saturation  level  at  hydrogen  inlet  is  that  the 
target  acid  concentration  inside  the  cell  is  around  90-95% 
H3PO4.  The  partial  pressure  of  water  over  this  acid  at  150  °C 
matches  that  of  the  water  vapor  pressure  at  60  °C.  In  order  to 
minimize  loss  of  moisture  at  hydrogen  side,  the  gas  is 
moistened  accordingly. 

Fig.  8  depicts  experimental  and  simulated  current  for  type 
A  grooves  at  different  oxygen  flowrates.  As  expected  with 


increase  in  flowrate  the  current  increases.  The  dotted  line 
shows  the  current  predicted  by  the  model  under  infinite 
oxygen  flow,  i.e.  if  a  single  dimension  model  without  the 
partial  pressure  drop  along  the  groove  is  considered. 

Fig.  9  shows  a  similar  current  vs  oxygen  flowrate  plot 
of  experimental  and  simulated  values  for  type  A  grooves 
at  a  lesser  overpotential  of  0.3  V  per  cell.  Here,  it  can  be 


Table  1 


Values  of  parameters 


Parameter  name 

Value 

Remark 

A* 

7iR2CTCT  (cm2)  of  projection  in  FZ-plane 

Literature 

aa 

3.5E— 04  (cm2  cm-3)  of  agglomerate 

Experimental 

b 

120  mV  per  decade 

Tuning 

Agg 

1.0E— 04  (cm2  s_1) 

Literature 

A 

0.2E— 03  (cm2  s-1) 

Literature 

Diffusion  layer 

400  pm  thick 

Experimental 

E 

0.92  V 

Experimental 

H 

5.5E— 08  (gmol  cm-3  atm) 

Literature 

lo 

1.0E— 06  (A  cm2) 

Literature 

n 

4 

Stoichiometry 

Agg 

1.0E— 04  (cm) 

SEM  picture 

Reaction  layer 

400  pm  thick 

Experimental 

T 

150  +  273  K 

Experimental 

W 

10  cm 

Experimental 

Ke 

Function  of  phosphoric 
acid  concentration 

Model 

Pa 

0.5/7iR2gg 

Literature 

2d 

0.5 

Experimental 

A 

0.5 

Experimental 

ac 

0.2  cm2 

Experimental 

ns 

13  grooves  per  manifold  (type  A) 

Experimental 
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seen  that  the  trend  is  similar  and  the  one-dimensional 
approach  would  have  predicted  a  higher  value  than  the 
actual. 

Figs.  10  and  1 1  show  the  moisture  taken  up  by  the  oxygen. 
It  can  be  seen  that  the  water  evaporated  is  much  higher  than 
the  water  generated  at  the  cathode  due  to  reaction.  The 
dotted  line  depicts  the  theoretical  water  generation  due  to 
reaction,  determined  from  the  total  experimental  current. 
It  can  be  also  seen  that  the  water  evaporation  from  the 
acid  layer  is  almost  independent  of  the  current  or  over¬ 
potential  and  depends  primarily  upon  the  gas  flow  rate. 


This  is  possibly  due  to  much  higher  water  uptake  of  oxygen 
than  water  generated  which  makes  the  outlet  gas  almost 
saturated. 

The  other  type  of  groove  pattern  on  which  the  validation 
runs  were  taken  was  “type  B”  grooves  as  shown  in  Fig.  5. 
The  trend  for  this  groove  pattern  remains  same  as  expected 
with  increase  in  current  with  flowrates.  Requirement  to 
explore  various  kinds  of  groove  patterns  and  their  efficacy 
comparison  with  the  help  of  the  steady  state  model  are 
discussed  in  the  Section  4.4.  Fig.  12  shows  current  vs 
flowrate  data  at  an  overpotential  of  0.3  V  for  type  B  grooves. 


Fig.  8.  Validation  data  for  steady  state  model — type  A  groove. 
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The  current  value  at  various  flowrates  of  oxygen  shows 
marginal  improvement. 

4.3.  Parameter  sensitivity  analysis 

The  complete  list  of  parameters  and  the  actual  values 
used  for  the  model  validation  is  referred  in  Table  1 .  Values 
of  some  of  these  parameters  like  diffusion  coefficients, 
vapor  pressure,  porosity  etc  were  taken  from  the  literature. 
Parameters  like  electrode  thickness,  catalyst  surface  area, 
agglomerate  diameter,  etc.  were  determined  experimentally. 
It  is  essential  to  determine  the  dependency  of  the  model 


output  on  the  accuracy  of  the  parameters.  Many  of  the 
parameters  used  in  this  model  cannot  be  determined  very 
accurately,  like  exchange  current  density,  active  area  of 
catalyst,  etc.  A  range  which  often  extends  in  decades  are 
possible  for  many  of  these  parameters.  Secondly,  for  tuning 
of  the  model,  to  represent  the  experimental  data,  parameter 
sensitivity  is  essential.  Apart  from  that,  estimating  output 
behavior  of  a  validated  model  provide  useful  insight  for  the 
actual  system  improvement. 

Of  the  parameters  used  in  the  model,  Tafel  slope  shows 
highest  sensitivity  towards  model  output.  In  order  to  test  this 
parameter  very  high  flowrates  of  oxygen  were  considered 


Fig.  10.  Moisture  loss  in  oxygen  side  at  0.4  V  ovepotential — type  A  groove. 
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which  eventually  makes  the  system  free  from  the  effects  of 
operating  conditions  like  gas  flow  rate  and  inlet  humidity. 
It  was  found  that  a  change  of  Tafel  slope  value  from  0.1 1  to 
0.12  V  per  decade  causes  the  output  current  density  to 
change  from  300  to  220  mA  per  sqcm  as  shown  in 
Fig.  13.  The  other  parameter  which  shows  strong  depen¬ 
dency  is  active  catalyst  area  per  unit  volume  of  agglomerate 
(af),  as  shown  in  Fig.  14.  The  effect  of  the  same  on  current 
output  is  almost  linear  as  expected. 

The  diffusion  resistance  inside  agglomerate  decreases 
with  agglomerate  size  resulting  in  increase  of  current  as 
shown  in  Fig.  15.  Dependency  of  current  on  agglomerate 


size  is  found  to  be  less  sensitive.  Oxygen  exchange  current 
density,  i.e.  the  rate  of  reaction  per  unit  area  of  catalyst 
surface  represented  in  the  form  of  current  at  equilibrium  is  a 
debated  one  as  it  strongly  depends  upon  the  nature  of  Pt 
crystal  structure,  presence  of  various  forms  of  platinum 
oxide,  etc.  [2,26].  Various  exchange  current  density  were 
tried  and  the  results  are  shown  in  Fig.  16. 

4.4.  Utility  of  the  model 

The  two-dimensional  steady  state  model  found  its  utility 
primarily  in  fuel  cell  component  design.  With  this  two- 


Fig.  12.  Validation  data  for  steady  state  model-  type  B  groove. 
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dimensional  model  one  can  analyse  the  various  sections 
of  an  electrode  and  can  study  the  variation  of  reaction  rates 
at  different  places.  This  helps  in  improvement  and  scale  up 
of  components  like  electrode  shape  factors,  reduction  of 
diffusional  resistance,  optimum  catalyst  layer  thickness,  etc. 

4.4.1.  Groove  design 

The  groove  pattern  on  graphite  manifold  plates  helps  in 
proper  distribution  of  gases  and  can  be  designed  for  better 
cell  performance.  For  example,  the  simple  type  A  groove 
causes  drop  in  stack  current  at  low  flowrates.  This  requires 


enhanced  flowrates  which  cause  loss  of  excess  oxygen.  To 
reduce  the  losses  one  can  split  the  gas  flow  into  a  narrower 
flow  width  and  increase  the  flow  path  length.  Type  B  groove 
pattern  is  a  step  towards  this  direction,  however,  a  large 
increase  in  path  length  may  cause  higher  hydraulic  pressure 
drop.  The  Figs.  9  and  12  show  marginal  improvement  of 
current  over  type  A  grooves. 

On  analysis  of  the  partial  pressure  contour  from  the 
model,  the  major  concentration  gradient  starts  at  the  inter¬ 
phase  of  catalyst  and  diffusion  layer.  To  overcome  this  a 
proprietary  “jump  over”  groove  design  is  developed  at 


Fig.  14.  Effect  of  active  catalyst  area  on  electrode  performance. 
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NMRL  and  analysed  with  this  model.  The  final  pattern 
yields  around  20%  higher  current  than  that  of  type  A  grooves 
at  low  flowrate. 

4.4.2.  Catalyst  layer  thickness 

As  the  catalyst  layer  thickness  increases  the  catalyst 
loading  and,  hence,  the  total  number  of  reaction  sites 
increase.  However,  with  increasing  thickness  the  resistance 
of  cation  to  flow  inside  the  catalyst  layer  increases.  The 
effect  of  catalyst  thickness  on  current  density  at  a  fixed 
overpotential  has  been  studied  using  the  model  which  shows 


that  a  catalyst  layer  of  around  400-500  pm  is  sufficient  for 
the  electrode  performance. 

4.4.3.  Humidity  control 

Control  of  cell  humidity  under  steady  state  operation  is 
generally  done  by  humidifying  the  hydrogen  gas,  slightly 
higher  than  the  equilibrium  water  vapor  pressure  of  the 
target  acid  concentration  at  the  fuel  cell  temperature. 
The  oxygen  is  generally  fed  dry.  This  creates  a  problem, 
as  the  amount  of  water  loss  from  acid  in  the  oxygen  stream  is 
quite  high  (refer  Fig.  10).  This  causes  problems  associated 


Fig.  16.  Effect  of  exchange  current  density  on  electrode  performance. 
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with  acid  drying.  With  the  help  of  the  model  the  water  loss 
can  be  predicted  and  humidification  of  hydrogen  gas  can  be 
done  accordingly. 

5.  Conclusions 

The  model  developed  helps  in  analyzing  the  distribution 
pattern  of  oxygen  partial  pressure  at  various  points  of  a 
composite  porous  gas  diffusion  cathode  of  a  PAFC.  Current 
generation  at  various  points  of  a  PAFC  cathode  depends 
upon  local  oxygen  concentration.  This  model  helps  in 
determining  the  area(s)  of  low  oxygen  concentration  and, 
thus,  is  helpful  for  designing  various  components  of  PAFC 
stacks.  The  model  has  been  extensively  used  in  NMRL  for 
designing  of  groove  pattern  for  the  oxygen  flow  field  in 
graphite  manifold  plates  and  the  optimum  catalyst  coating 
thickness. 
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